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Complex systems displaying recurrent spike patterns are ubiquitous in nature. Understanding the 
organization of these patterns is a challenging task. Here we study experimentally the spiking output of a 
semiconductor laser with feedback. By using symbolic analysis we unveil a nontrivial organization of 
patterns, revealing serial spike correlations. The probabilities of the patterns display a well-defined, 
hierarchical and clustered structure that can be understood in terms of a delayed model. Most importantly, 
we identify a minimal model, a modified circle map, which displays the same symbolic organization. The 
validity of this minimal model is confirmed by analyzing the output of the forced laser. Since the circle map 
describes many dynamical systems, including neurons and cardiac cells, our results suggest that similar 
correlations and hierarchies of patterns can be found in other systems. Our findings also pave the way for 
optical neurons that could provide a controllable set up to mimic neuronal activity. 

Nature presents many fascinating complex systems in which distinguishing signatures of determinism in 
their high-dimensional dynamics is extremely challenging 1 4 , not only because of the presence of noise, 
but also, because of lack of information about the system: one can measure only one or few relevant 
variables, and with a limited spatial and/or temporal resolution. A successful approach for studying such systems 
is by focusing on an event-level description of their dynamics, considering, for example, intervals between events. 
Examples of this approach include neuronal inter-spike intervals, heart beat-to-beat intervals, earthquake waiting 
times, intervals between peak communication activity in social networks, etc. 510 . For the analysis of these events 
the symbolic method known as ordinal analysis 1114 considers the relative order in which the events occur. 

Here we use ordinal analysis to study the spiking output of a stochastic optical system, consisting of a 
semiconductor laser with feedback from an external reflector. Close to the lasing threshold, moderate feedback 
levels induce apparently random spikes in the laser output intensity, which become more frequent as the laser 
pump current is increased (see Figs, la and lb). In spite of the fact that the spiking dynamics (referred to as low- 
frequency-fluctuations, LFFs) has attracted a lot of attention in the last decades 15 , several features of the dynamics 
remain poorly understood, in particular, the presence of correlations in the spike sequence, and how these 
correlations change as a control parameter, such as the laser pump current, is varied. The well-known Lang 
and Kobayashi (LK) model 16 predicts that, for typical parameters, the LFF is a transient dynamics sustained by 
spontaneous emission noise, and the duration of the deterministic transient increases with the pump current 17,18 . 
Therefore, the presence of correlations in the spike sequence is the fingerprint of the underlying attractor. The 
main goals of our work are to characterize the symbolic dynamics underlying the spike sequence in order to infer 
spike correlations, to test whether simulations of the stochastic LK model also display similar correlations, and to 
investigate if a similar behavior can be observed in other natural systems; specifically, we aim to answer the 
question whether in sequences of optical spikes there are serial correlations similar to those present in spike 
sequences of biological neurons. 

With these aims we analyze experimentally recorded sequences of optical spikes and unveil a hierarchical and 
clustered organization of ordinal patterns (in the following referred to as words), with a well defined structure of 
the probabilities of occurrence. Simulations of the LK model are found to be in good agreement with the 
observations. To the best of our knowledge, ours is the first demonstration of an organized symbolic structure 
underlying the LFF dynamics, revealing serial correlations among several consecutive intensity dropouts. 

Most importantly, we identify a minimal iterative model, a modified circle map (proposed in 19 to represent 
spike correlations in sensory neurons) that, despite its simplicity, displays a symbolic dynamics with the same 
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Figure 1 | Spiking laser output, experimental setup, and word probabilities, (a, b), Experimental time trace of the laser intensity when the pump current 
is 26.5 mA and 27.3 mA respectively. The red dots indicate the spike times and a few ordinal patterns (words) are indicated as examples, (c), Experimental 
setup, where LD stands for laser diode, BS for beamsplitter, PD for photodetector, OSC for oscilloscope and GRAT for grating (see Methods), 
(d), Probabilities of the six words (schematically indicated in the right panel) vs the laser pump current. Two clusters of words can be observed, with 
similar probabilities: 021-102 (black-green) and 120-201 (magenta-cyan). A crossover in the hierarchical organization of the words occurs at about 
26.6 mA: at lower current values the word 012 (blue) is the most probable one, while at higher currents values, the word 210 (red) is the most probable 
one. The error bars are estimated with a binomial test and the gray region indicates probability values consistent with the null hypothesis that there are no 
correlations in the spike sequence, and therefore, the six words are equally probable 19 . 



hierarchical and clustered organization of patterns as the experi- 
mental data. In order to confirm that this minimal model indeed 
represents the symbolic dynamics underlying the sequence of optical 
spikes, we consider the spiking output of the laser with periodic 
forcing (via a direct modulation of the pump current). We dem- 
onstrate that the symbolic dynamics of the forced laser is also in 
good qualitative agreement with the symbolic dynamics of the modi- 
fied circle map. 

Results 

The experimental setup (schematically shown in Fig. lc and de- 
scribed in Methods) uses a commercial, inexpensive semiconductor 
laser. We recorded long sequences of optical spikes and analyzed the 
inter-spike intervals, {AT,- = f I+1 — t,-}, with f, being the time when a 
spike occurs. The data contains 70,000-300,000 spikes, at low and 
high pump current respectively. In the Supplementary Information 
the threshold criterium used to define the spike times is discussed, 
and the underlying symbolic organization is shown to be robust to 
threshold variations. 

We construct the ordinal patterns (words) with D = 3 consecutive 
intervals (four consecutive optical spikes) as described in Methods. 
In Figs, laand lbafewwords are indicated as examples. Results for D 
= 2 (and D = 4) are presented in the Supplementary Information. 
Figure Id displays the probabilities of the six possible words (indi- 
cated in the right panel of Fig. Id), computed from time-series 
recorded within a wide range of pump currents. The symbolic ana- 
lysis reveals that serial correlations among four consecutive spikes 
are present in the spike sequence: if there are no correlations (null 
hypothesis, represented as gray region in Fig. Id), the six words are 
equally probable; however, in Fig. Id one can observe that the prob- 
abilities are clearly not consistent with the null hypothesis. 

Moreover, one can observe a hierarchy in the probability values, 
which presents a crossover at about 26.6 mA: for higher pump 



currents the most probable word is 210 (corresponding to three 
consecutively decreasing inter-spike intervals), while for lower pump 
currents, the most probable word is 012 (corresponding to three 
consecutively increasing intervals). 

Figure Id also reveals a clustered organization of the probabilities: 
words 021 and 102, on the one hand, and words 120 and 201, on the 
other hand, occur with similar probability. The probabilities of these 
two pairs of words present the same evolution when the pump cur- 
rent is varied. The robustness of these findings was confirmed by 
using different lasers and feedback conditions (see Supplementary 
Information), where the same clustered hierarchical structure was 
found; in fact it can also be observed in Fig. 2b of. 

To further check the robustness of our observations, we performed 
simulations of the LK model 16 (see Methods). The model predicts a 
fast pulsing behavior, in the picosecond time-scale that, after being 
properly filtered out (to account for the bandwidth of experimental 
measurements), gives a sequence of spikes in good agreement with 
the observations. 

Figures 2a and 2b show the numerically calculated word probabil- 
ities vs the laser pump current parameter of the model, ji, for two sets 
of parameters. The simulations agree qualitatively well with the 
experimental findings: we observe the same hierarchical and clus- 
tered organization of patterns, as well as the same crossover. 
Figures 2a and 2b also point to the robustness of the results, as the 
two clusters are present in simulations with different feedback 
strengths. 

A relevant question is whether the hierarchical and clustered 
organization of symbolic patterns uncovered here also occurs in 
other natural systems. If this is the case, there should be a minimal 
model that also displays such organization of ordinal patterns. We 
analyzed the symbolic dynamics of several iterative maps and found 
one, surprisingly simple, a modified circle map (see Methods), that 
reproduces the hierarchy and clustered organization of the word 
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Figure 2 | Probabilities of the ordinal patterns calculated numerically, (a, b), The probabilities are computed from simulations of the LK model; the 
feedback strength is r\ = 10 ns~' and 20 ns~' respectively; other model parameters are as indicated in Methods, (c), Probabilities computed from the 
minimal model, Eq. (3). The parameters are: p = 0.23, K = 0.04 and ft c = 0.002. Both, in the LK model and in the minimal model, the hierarchy in the 
word probabilities, the clustering and the crossover are the same as in the experimental data, Fig. Id. 



probabilities. Figure 2c shows the probabilities obtained from the 
modified circle map vs the map parameter, a c . We note that the effect 
of varying a c corresponds, qualitatively, to varying the laser pump 
current: a similar behavior as in Figs. Id, 2a and 2b can be observed. 
There is the same hierarchy in the probabilities, which undergoes the 
same crossover. The existence of the two clusters is also remarkable, 
formed by the same words as in the sequence of optical spikes. 

To further demonstrate that the modified circle map is a minimal 
model for the spiking laser output, we consider the laser response to a 
periodic forcing. Figures 3a, 3b display typical time traces of the laser 
output with direct pump current modulation, and Fig. 3c displays the 
experimental setup (see Methods). In Fig. 3a, the effect of weak 
forcing can be clearly observed and, as the forcing increases, the 
spikes synchronize with the forcing and the inter-spike intervals 
become multiples of the modulation period, Fig. 3b 21,22 . 

Figure 4a displays the word probabilities computed from the 
experimental data, as a function of the modulation amplitude. We 
observe that, despite the external forcing acting on the system, the 
symbolic dynamics has the same two clusters as in the unperturbed 
case (021-102 and 120-201). Also a clear hierarchy is manifested in 
the probabilities, which present a crossover too, at about 2% of 
modulation amplitude. For strong modulation the hierarchy of pat- 
terns is modified, with words 012 and 210 now being the less prob- 
able ones. This can be understood by considering that the 
modulation imposes a "rhythm" in the spike sequence and consecu- 
tively decreasing or increasing inter-spike intervals are less probable 
as they imply that the synchrony with the external rhythm is lost, 
either because the spikes are increasingly close (AT t > AT 1+1 > 
AT j+2 representing word 210), or because the spikes are increasingly 
separated (AT, < AT i+1 < AT 1+2 , word 012). Simulations of the LK 
model including a sinusoidal modulation of the pump current para- 
meter were found to be in good agreement with the observations, and 
a detailed comparison model-experiments will be reported 
elsewhere. 

The symbolic behavior is in qualitative good agreement with the 
modified circle map, when the parameter K, that represents the 
forcing amplitude, is varied (Fig. 4b). However, the agreement is 
not as good as in the unperturbed situation and we interpret that 
this could be due to the fact that in the experiment, the current 
modulation introduces an additional source of noise, which probably 
depends on the modulation amplitude. In order to test this inter- 
pretation we considered an additional noise term in the map equa- 
tion (see Methods) and, as can be seen in Fig. 4c, this term mainly 
affects the word probabilities at low modulation amplitudes (before 
the crossover), improving the agreement with the experimental 
probabilities (Fig. 4a). As could be expected, the symbolic behavior 
depends not only on the modulation amplitude, but also, on the 



modulation frequency, and a detailed analysis is in progress and will 
be reported elsewhere. 

Discussion 

The symbolic behavior of the modified circle map demonstrates that, 
with and without periodic forcing, it is an adequate minimal model 
for representing qualitatively the symbolic dynamics underlying the 
sequence of optical spikes. Considering the high-dimensionality of 
the laser dynamics (induced by the feedback delay time) it is unex- 
pected and remarkable that a simple iterative model with only one 
variable can mimic the extended correlations present in the sequence 
of spikes. 

Both, in the experiments and in the LK model simulations, we 
found the presence of two clusters of word probabilities (021-102, 
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Figure 3 | Laser spiking output under external periodic forcing. Time 
trace of the laser intensity when the modulation amplitude is (a), 31.2 mV 
(1.6% of the DC current); (b), 70.2 mV (3.6% of the DC current); other 
parameters are as indicated in Methods. The effect of the forcing is clearly 
visible in the laser intensity, (c), Schematic representation of the 
experimental setup, where the forcing is included as a direct modulation of 
the laser pump current. LD stands for laser diode, WG for waveform 
generator, BS for beamsplitter, PD for photodetector, OSC for oscilloscope 
and M for mirror (see Methods). 
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Figure 4 | Influence of periodic forcing in the organization of the symbolic patterns, (a), Probabilities computed from the experimental data. The DC 
value of laser pump current is 39 mA and the modulation amplitude is increased from 0 mVto 78 mV (0% to 4% of the DC current, other experimental 
parameters are as described in Methods). The two word clusters are clearly seen, which are the same as without modulation. A crossover in the word 
hierarchy can also be observed at about 1.8%. (b), Probabilities computed from the modified circle map, when varying the parameter iCwith constant 
noise strength. Here again, the hierarchical and clustered structure of the probabilities is in good qualitative agreement with that in the experimental data. 
The parameters are p = —0.23, a c = 0.2 and ji c = 0.002. (c), Probabilities computed with the modified circle map considering a noise strength that 
decreases with the modulation amplitude. For weak modulation the agreement with the experimental probabilities is further improved. 



120-201). We also found that 012 and 210 are the most probable 
words, at low and high current respectively, revealing that serial 
correlations, among four consecutive optical spikes, vary with the 
laser pump current. We show in the Supplementary Information that 
these correlations are not detected via the standard correlation ana- 
lysis. The good qualitative agreement experiment-simulations indi- 
cate that this organization of symbolic patterns is a fingerprint of the 
topology of the attractor underlying the LFF dynamics. 

As a future study it would be interesting to consider the phe- 
nomenological model proposed by Mindlin and coworkers 23 , which 
gives a good description of the probability distribution function of 
the inter-spike intervals. In the symbolic dynamics of this model the 
most probable words are 012 and 2 1 0 24 , and it would be interesting to 
explore the parameter space to search for the two clusters of words 
uncovered here. 

To summarize, we demonstrated the existence of a hierarchical 
organization of patterns in the symbolic dynamics of semiconductor 
lasers, in the spiking regime induced by delayed optical feedback. 
This organization displays a clustered hierarchical structure which 
had not been previously noticed, despite the great attention that the 
laser spiking output has attracted. Simulations of the LK model were 
found to be in good qualitative agreement with the experimental 
observations, providing in this way a validation of the LK model in 
unprecedent long time-scales. We also found a minimal model, a 
modified circle map, that describes the hierarchy and clustered sym- 
bolic structure. The validity of the minimal model was further 
demonstrated by analyzing the symbolic dynamics of the laser with 
periodic forcing via direct current modulation. 

Since the modified circle map is a generic model of forced oscilla- 
tors, we envision that the underlying symbolic structure found here 
could be observed in many other dynamical systems. In particular, 
our results could be relevant for improving our understanding of the 
firing patterns of biological neurons. Establishing a direct connection 
between these different dynamical systems can offer new perspec- 
tives, both, in photonics and in neuroscience. Since the modified 
circle map has been proposed as a minimal model of electroreceptors 
of paddlefish 18 , our results suggest that optical neurons inspired in 
biological ones, displaying similar temporal correlations in their fir- 
ing patterns, could be built using semiconductor lasers, thus provid- 
ing a novel, inexpensive and controllable experimental set up for 
mimicking neuronal activity. In particular, the optical setup could 
be used to analyze the role of external forcing in the spike sequence, 
providing, for example, new insight in the role of gamma oscillations 
in the brain, which have been shown to serve to concentrate neuronal 



discharges to particular phases of the gamma cycle 25 . On the other 
hand, because semiconductor lasers are nowadays the coherent light 
sources used in telecommunications and data storage, and because 
there is a wide range of spiking regimes with experimentally control- 
lable parameters 26 , we envision that our finding could also contribute 
to the development of novel applications of semiconductor lasers in 
neuro-inspired optical computing devices 27 " 29 . 

Methods 

Experimental setup. The experiments were performed with a 675 nm Al-GalnP 
semiconductor laser (AlGalnP Sanyo DL-2038-023) with optical feedback from a 
diffraction grating (see Fig. lc). The external cavity length was 70 cm, giving a 
feedback delay time of 4.7 ns. A beam-splitter sends 50% of the light to a 1 GHz 
oscilloscope (Agilent Infmiium 9000). The laser temperature and pump current were 
controlled to an accuracy of 0.01 C and 0.01 mA respectively with a ITC502 Thorlabs 
laser diode combinator controller. The operating temperature was 18 C and the laser 
pump current was varied in steps of 0.1 mA, from 26.3 mA to 27.3 mA. At 18 C the 
threshold current of the solitary laser is l th = 27.8 mA, and the feedback-induced 
threshold reduction is 6.5%. Time series of 32 ms were recorded (3.2 X 10 7 data 
points), containing 70,000-300,000 spike events, for low and high pump current 
respectively. 

To study the dynamics under periodic forcing, a Hitachi HL6714G diode laser was 
used, lasing at 650 nm, with solitary threshold of /(/, = 38 mA. The pump current was 
set to 39 mA, and optical feedback resulted in a threshold reduction of 8%. The 
external cavity length was of 70 cm, giving a time delay of 4.7 ns. The pump current 
was modulated through a bias-tee in the laser mounting and with an 80 MHz 
waveform generator (Agilent 33250A). The modulation frequency was 17 MHz and 
the modulation amplitude was varied in the range 0% to 4% (78 mV), relative to the 
DC value of 39 mA (1.037^). Time series of 20 ms were recorded, which contained 
between 70,000 and 200,000 spike events, for low and high modulation amplitude 
respectively. 

Ordinal time-series analysis. The time series of inter-spike intervals, {AT;}, were 
analyzed using the ordinal method 11 , which takes into account the order in which the 
values appear in the time series. As this method is well known, we only describe it 
briefly. 

A time series {x(t)} of length N is divided in N — D vectors of length D, and each 
element of a vector is replaced by a number from 0 to D — 1, in accordance with the 
relative length of the element in the ordered vector (0 corresponding to the shortest 
and D — 1 to the longest value in each vector). In this way, each vector has associated 
an "ordinal pattern" (or word) composed of D symbols. For example, with D = 3, the 
vector {x(t), x(t + 1), x(t + 2)} = {5, 1, 10} gives the word 102, because x(t + 1) < x(t) 
< x(t + 2). Then, by counting the frequency of occurrence of the different words, 
their probabilities are computed. To perform an statistically significant analysis, the 
number of data points in the time series should be N ^> D\. This condition is fulfilled 
as we perform the analysis with D = 2, 3, and 4, and in our datasets N is between 
70,000 and 300,000. 

In all the figures the probabilities of the words have been computed considering a 
rolling-window. Specifically, for D = 3 (D = 2) the time-series of inter-spike intervals, 
{ATi, AT 2 . . ., AT N }, was transformed into 3 (2) symbolic sequences by considering 3 
(2) different initial ATvalues for encoding the words. For example, for D = 3, in one 
sequence the words are formed as (AT 1 ,AX 2 ,AT 3 ), (Ar 4 ,AX5,AT 6 ), etc. In the second 



SCIENTIFIC REPORTS | 4:4696 | DOI: 1 0. 1 038/srep04696 



4 



sequence as (AT 2 ->ATj, > AT 4 ), (AT 5 ,AT 6 ,AT 7 ), etc. And in the third sequence as 
(AX 3 ,AT 4 ,AT'5), (AT 6 ,AT 7i AT 8 ), etc. From each of these sequences, the probabilities 
of the words were computed, an then, they were averaged before plotting the results. 
This is fully equivalent to consider overlapping patterns (i.e., a single symbolic 
sequence formed as (AT U AT 2 , AT 3 ), (AT 2 , AT 3 , AT 4 ), etc.) 

Compared with other symbolic methods, the ordinal transformation has the 
advantage that the symbolic sequence keeps the information about the correlations 
present in the time-series, but it has the drawback that the information about the 
actual values (i.e., the precise duration of the inter-spike intervals) is lost. However, 
this drawback is in fact an advantage in our case, as it allows using a simple thresh- 
olding method to determine the spike times. Because the spikes have fluctuating 
amplitude (particularly at low pump currents, see Fig. la), a simple threshold might 
give rise to timing errors; however, with ordinal analysis the precise values of the 
intervals are disregarded (as only the relative order of consecutive intervals is con- 
sidered); therefore, with the ordinal transformation there is no need of using methods 
more sophisticated than thresholding to precisely determine the timing of the events. 
While the threshold chosen determines the number of spikes detected, and thus, the 
symbolic sequence of ordinal patterns obtained, in the Supplementary Information 
we also demonstrate that our findings are robust, with the hierarchical and clustered 
structure of word probabilities persisting in a wide-range of threshold values. 

The LK model. The Lang and Kobayashi (LK) model 15 consists of two coupled delay- 
differential rate equations governing the evolution of the slowly varying complex 
electric field amplitude, E, and the carrier density, JV. The model equations are: 

dE 1 / 

dt = i7„ (1 +a)(G ~ l ) E +nE{t-^)e- w + y/a^f, (l) 

__(„_ W _ G | £ | 2) , (2) 

where x p and t n are the photon and carrier lifetimes respectively, a is the linewidth 
enhancement factor, G is the optical gain, G — N j ( 1 + e|£| ) (with € being a 
saturation coefficient), fj, is the pump current parameter, ;/ is the feedback strength, t 
is the feedback delay time, co 0 z is the feedback phase, fi sp is the noise strength, 
representing spontaneous emission, and £ is a Gaussian distribution with zero mean 
and unit variance. 

Typical values of the parameters were used in the simulations: k — 300 ns" 1 , ^ — 
1 ns 1 , € — 0.01, t — 5 ns,and/? s p= 10~ 4 ns -1 ; in Figs. 2a and 2b the pump current, ji, 
is the control parameter and the simulations were performed with two sets of para- 
meters which were found to give good agreement with the experimental data (in 
Fig. 2a, i] — 10 ns -1 , x = 4; in 2b, i] = 20 ns -1 , a — 4.5). Because the LK model is a 
simple model (it takes into account only one reflection in the external cavity, it 
neglects multi-mode emission, spatial effects, and thermal effects; in particular the 
shift of the emission wavelength with increasing current due to Joule heating), only a 
qualitative agreement could be expected. 

Within the framework of the LK model it has been shown that the optical spikes are 
often a transient dynamics, even when noise is included in the simulations, and the 
duration of the transient increases with the pump current parameter 17 . In order to 
compute long time series of inter-dropout intervals, when the trajectory found a 
stable fixed point the simulation was re-started with new random initial conditions. 
Simulations of 2 milliseconds were performed, obtaining time series containing 
between 12,000 and 48,000 events, for low and high pump currents, respectively. 

Minimal model. The modified circle map proposed by Neiman and Russell 19 to 
describe serial correlations in electroreceptors of paddlefish was found to display the 
same symbolic organization as that found of the spiking output of the laser with 
feedback, with and without periodic forcing. The map equation is: 

^(< + l)=(*(i) + p+^[sin(2^(i))+« c sta(47^(i) )] + &{(*) (3) 
In 

The parameter p is proportional to the ratio between the period of the forcing and the 
natural period of the oscillator, K is proportional to the forcing amplitude, a c 
represents the strength of a second harmonic component, and fi c represents the 
strength of stochastic fluctuations (q being a Gaussian white noise). Considering that 
the values <j)(i) can represent the spike times, we analyzed the time series of phase 
increments, <j)(i + 1) — <f)(i), in analogy with the inter-spike intervals, AT,. As Figs. 2c 
and 4c show, varying the parameters c/. c and K correspond to varying the laser pump 
current and the modulation amplitude, respectively. To take into account the 
additional noise introduced by the external forcing, a noise term inversely 
proportional to the modulation amplitude [/? e (l + l/K)c(i)] is considered in Fig. 4d. 
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